function [savings_delta] = generate_delta_CI_table(covariate,wave,winter_prevuse,SMC,fbase,scaled_ATT,s_now_batch,rule,tcost)
% This function summarizes point estimate and confidence interval of 
% the difference between EWM savings and RCT benchmark 

% Inputs:
% (1) covariate: indicates covariate for analysis: income, size, vintage,
% minimum of baseline consumption, maximum of baseline consumption, or
% standard deviation of consumption
% (2) wave: indicates whether the output is wave-specific (=3,6,or 7) or
% for the pooled sample (=0) 
% (3) winter_prevuse: indicates whether baseline consumption is calculated 
% as the mean of consumption in winter months (Jan and Feb) or as the mean
% of specified pre-treatment periods
% (4) SMC: =1 use social marginal cost; =0 use retail electricity price
% (5) fbase: string used to indicate the propensity score and baseline
% months specification for output table 
% (6) scaled_ATT: scaled att as a benchmark 
% (7) s_now_batch: timestamp used to identify EWM bootstraps 
% (9) rule: treatment rule 
% (10) tcost: private marginal cost of implementing the program per household
%  >0 represents cost savings; =0 represents kWh reduction 

% Outputs:
% (1) savings_delta: summary table of difference between EWM savings and 
% RCT benchmark

%% obtain point estimates 
switch winter_prevuse 
    case 1
        s_winter = '_winter';
    case 0
        s_winter = '';
end
if (tcost)>0 
    if (SMC)
        s_cost = 'SMC';
    else
        s_cost = 'PMC';
    end
elseif (tcost)==0
    s_cost='kwh';
end

if rule=="quadrant"
    filename_coefs=sprintf('coef_quadrant_%s_%s%s_wave%1.0f.mat',covariate,s_cost,s_winter,wave);
elseif rule=="cubic"
    filename_coefs=sprintf('coef_cubic_%s_%s%s_wave%1.0f.mat',covariate,s_cost,s_winter,wave);
elseif rule=="onedim"
    filename_coefs=sprintf('coef_onedim_%s_baseline_%s%s_wave%1.0f.mat',covariate,s_cost,s_winter,wave);
end
    file = sprintf('%s_%s',fbase,filename_coefs);
S=load(file);

%extract content of S_quadrant
c_fieldnames = fieldnames(S);
for ifield = 1:length(c_fieldnames)
    field_ = c_fieldnames{ifield};
    eval(sprintf('%s=S.%s;',field_,field_))
end

% obtain scaled att
scaled_att_PE=scaled_ATT;
%% obtain ewm point estimate (PE) and confidence interval (CI)
% generate suffix for bootstrap folder 
if (tcost)
    if (SMC)
        s_speci = 'smc';
    else
        s_speci = 'pmc';
    end
else
    s_speci='kwh';
end
if rule=="quadrant"
    fpath_BS = ['.\intermediate_data\EWM_results\BS\quadrant\',...
     s_now_batch,'_',covariate,'_',s_speci,'_','deltaCI','\'];
elseif rule=="cubic"
     fpath_BS = ['.\intermediate_data\EWM_results\BS\cubic\',...
     s_now_batch,'_',covariate,'_',s_speci,'_','deltaCI','\'];
elseif rule=="onedim"
    fpath_BS = ['.\intermediate_data\EWM_results\BS\quadrant\',...
        s_now_batch,'_','univariate','_',s_speci,'\'];
end
if ~exist(fpath_BS,'dir')
   mkdir(fpath_BS); 
end

[bs_result] = pool_BS_delta_CI(s_now_batch,rule,fpath_BS);

%% hypothesis testing
% Calculate PE of V(G_hat) for savings and CI for savings 
if rule=="quadrant" || rule=="onedim"
    percent = sum(nw.*in_Ghat)/n;
    savings = sum(gu.*in_Ghat)/n;
elseif rule=="cubic"
    percent = mean(in_Ghat);
    savings = nanmean(g.*in_Ghat)*Yscale;
end

% calculate delta difference between ewm point estimates and scaled att
delta_ewm_rct = savings - scaled_att_PE; %PE of difference on the original sample

% get bootstrap V_Ghat and scaled_att
V_Ghat = bs_result.W_bs_Ghat; 
scaled_att = bs_result.scaled_ATT_bs;

% calculate CI
se_delta=sqrt(var(V_Ghat-scaled_att));
ci_delta_lb = delta_ewm_rct - 1.96*se_delta;
ci_delta_ub = delta_ewm_rct + 1.96*se_delta;

%% Output for table 
savings_delta=nan(1,3);
savings_delta=array2table(savings_delta,'VariableNames',{'Covariates','Share\ treated','Difference\ in\ savings\ between\ EWM\ RCT'}); 
format short g
savings_delta.(1) = {covariate};
savings_delta.(2) =round(percent*100,1);
savings_delta.(3) =round(delta_ewm_rct,2);

ci_delta_lb=round(ci_delta_lb,2);
ci_delta_ub=round(ci_delta_ub,2);

ci_table=nan(1,3);
ci_table=array2table(ci_table,'VariableNames',{'Covariates','Share\ treated','Difference\ in\ savings\ between\ EWM\ RCT'}); 
ci_delta_lb_string=string(ci_delta_lb);
ci_delta_ub_string=string(ci_delta_ub);
ci_delta_for_table=strcat('(',ci_delta_lb_string,',',ci_delta_ub_string,')');
ci_table.(1)="";
ci_table.(2)="";
ci_table.(3)=ci_delta_for_table;

savings_delta=[savings_delta;ci_table];

end

